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A  NONLINEAR  ELLIPTIC  PROBLEM  RELATED  TO  FLOWING 
GRANULAR  MATERIALS 

PIERRE  A.  GREMAUD* 


Abstract.  Similarity  solutions  for  the  flow  of  granular  materials  are  constructed.  Unlike  pre¬ 
vious  work,  the  present  approach  can  be  applied  to  non-axisymmetric  containers.  The  steady  state 
equations  are  reduced  to  a  nonlinear  Helmholtz  equation  on  a  subdomain  of  the  sphere.  Exis¬ 
tence  and  local  uniqueness  of  the  solutions  are  established.  A  spectral  numerical  method  using 
Fourier/ Chebyshev-Gauss-Radau  collocation  for  discretization  and  Newton- GMRES  as  solver  is  pro¬ 
posed  and  implemented.  Corresponding  numerical  experiments  are  discussed. 
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1.  Introduction.  This  paper  is  about  the  determination  of  the  flow  of  granular 
material  under  gravity  in  hoppers  of  simple  geometry.  This  problem  is  common  to 
many  industrial  processes.  In  spite  of  its  apparent  simplicity,  it  presents  a  formidable 
array  of  difficulties.  In  practice,  the  withdrawal  of  stored  materials  from  hoppers  and 
bins  is  well  known  to  be  problematic.  No  flow,  segregation,  flooding  (uncontrolled 
flow)  and  structural  failures  are  often  encountered  [11].  Improved  design  criteria  are 
sought  through  a  better  understanding  of  such  flows.  All  the  previous  contributions 
we  are  aware  of  [4],  [5],  [7]-[9],  [13]-[17],  [20],  to  cite  but  a  few,  deal  only  with 
axisymmetric  containers.  This  significant  restriction  is  removed  here. 

Apart  from  the  geometry  of  the  hopper,  two  important  factors  are  the  internal 
friction  of  the  material  and  the  friction  between  wall  and  material.  Those  parameters 
are  further  discussed  in  Section  2.  The  flows  themselves  can  reach  from  established 
mass  flows  where  the  material  moves  and  deforms  everywhere  in  a  “smooth  way”  to 
highly  time-dependent  funnel  flows  where  motion  only  takes  place  in  the  central  part 
of  the  silo  [14].  In  the  first  case,  steady  state  models  are  often  appropriate.  Indeed, 
current  design  criteria  are  based  on  mass  flow  calculations  through  the  resolution  of 
steady  state  models,  a  view  point  that  is  also  taken  in  this  paper.  Throughout,  the 
material  is  assumed  to  be  an  incompressible,  perfectly  plastic,  cohesionless  Coulomb 
powder  with  a  yield  surface  of  von  Mises  type.  Further,  the  eigenvectors  of  the  strain 
rate  and  stress  tensors  are  assumed  to  be  parallel.  Those  assumptions  are  standard  in 
this  field,  although  the  alignment  condition  is  somewhat  controversial.  Those  issues 
are  commented  on  in  Section  2.  Apart  from  the  above  constitutive  relations,  the 
equations  derive  from  the  basic  principles  of  Continuum  Mechanics:  conservation 
of  mass  and  momentum.  Conservation  of  energy  does  not  enter  the  problem  as  heat 
losses  through  friction  influence  the  temperature  which  is  not  related  in  any  significant 
way  to  the  variables  under  study. 

Much  remains  to  be  done  about  the  analysis  of  the  above  problems  and  as  a  result, 
the  design  of  mass  flow  hoppers  still  greatly  relies  on  an  observation  due  to  Jenike 
[8],  [9]  in  the  late  19507s.  It  was  noticed  that  for  three  dimensional  conical  hoppers  of 
circular  cross  section  as  well  as  for  two-dimensional  wedge  shaped  hoppers  the  above 
equations  admit  similarity  solutions .  Those  radial  solutions  correspond  to  particle 
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paths  that  are  radial  lines.  The  similarity  is  reflected  in  scalings  of  the  stresses  and 
velocities  with  respect  to  the  radial  distance  r  from  the  vertex  of  the  hopper.  In  those 
simple  geometries,  the  equations  for  the  radial  fields  reduce  to  systems  of  ODEs  that 
can  be  easily  numerically  integrated.  The  stability  of  the  solutions  to  perturbations, 
change  of  the  geometry  (opening  angle)  and/or  physical  properties  of  the  hoppers  were 
studied  in  [4],  [13],  [15]  and  [16]  through  numerical  resolutions  of  the  full  PDE  system 
(in  [16],  the  material  is  assumed  to  be  compressible)  and  in  [5]  and  [15]  through  linear 
stability  studies.  Experimental  evidences  [14]  confirm  the  important  role  played  by 
the  radial  fields  in  practice. 

The  purpose  of  this  paper  is  to  show  that  Jenike’s  construction  of  similarity  so¬ 
lutions  can  be  generalized  to  general  mathematical  cones,  i.e.,  pyramidal  domains  of 
arbitrary  cross  section.  This  is  significant  as  previous  works  in  this  area  deal  exclu¬ 
sively  with  axisymmetric  containers,  even  though  those  are  the  exception  rather  than 
the  rule  in  practice  (see  e.g.  [18]  for  remarks  on  the  influence  of  the  hopper  geometry 
on  the  flowing  properties).  The  loss  of  axisymmetry  considerably  complicates  the 
structure  of  the  stress  tensor  and  the  ensuing  equations.  As  is  the  case  for  Jenike’s 
radial  solutions,  the  radial  structure  leads  to  significant  simplifications.  However,  in¬ 
stead  of  ODEs,  the  radial  stress  field  must  here  be  obtained  through  the  resolution  of 
a  nonlinear  Helmholtz  equation  in  a  subdomain  of  the  sphere. 

The  paper  is  organized  as  follows.  The  model,  geometry  and  physical  assumptions 
are  discussed  in  Section  2.  An  existence  and  local  uniqueness  result  is  stated  and 
proved  in  Section  3.  In  Section  4,  a  pseudospectral  numerical  method  is  proposed. 
It  uses  Fourier  collocation  in  longitude  and  Chebyshev-Gauss-Radau  collocation  in 
latitude,  to  account  for  the  boundary  conditions  at  hand.  Section  5  is  devoted  the 
description  of  several  numerical  experiments.  Conclusions  are  offered  in  Section  6. 

2.  The  model.  The  physical  quantities  and  corresponding  equilibrium  equations 
are  expressed  in  spherical  polar  coordinates,  with  the  origin  corresponding  to  the 
vertex  of  the  hopper.  For  non-axisymmetric  domains,  coordinate  systems  that  are 
better  suited  to  the  geometry  at  hand  can  usually  be  constructed.  However,  such 
systems  are  typically  not  orthogonal,  complicating  greatly  the  structure  of  the  basic 
equations  of  Continuum  Mechanics  [1].  To  simplify  the  numerics,  such  alternate 
coordinate  systems  are  introduced  in  Section  4  through  a  change  of  variables,  but  the 
individual  components  of  the  stress  tensor,  velocity,  etc...,  are  still  measured  in  terms 
in  the  original  spherical  coordinates. 

The  strain  rate  tensor  V  -  -$(Vv  +  Vur),  v  being  the  velocity,  and  the  stress 
tensor  T  take  respectively  the  form  (see  e.g.[19],  p.184)1 


V  = 


- drvr  -  \  (\9qVt  “  t  +  dTve)  -  -r  +  drVf) 

-$(vr  +  deve)  -jp  {dev<t>  -  cot <9^  +  ^b<V*) 
-±(^r  +  cot0U0  +  ^d0^) 


T  = 


Trr  Tr$ 
Tee 


The  equations  of  equilibrium  are 

(2.1)  V  •  T  =  pg, 


i\Ve  omit  to  write  the  lower  triangular  part  of  symmetric  tensors. 
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where  g  is  the  acceleration  vector  due  gravity,  or  equivalently 

drTrr  4-  — 4-  -deTro  4-  -(2 Trr  -  TH  -  Tee  4-  Ttq  cot  Q)  =  -pgcosd, 
rsrnfl  r  r 

drTr(h  4 - : — -c^T^  4 — OqT^q  4 — (3Tr^  4-  2T^  cot#)  =  0, 

v  rsmOY  r  r 

drTre  +  — +  I  (3Tre  +  ( Tgg  -  T^)  cot  0)  =  pg sin  1 9, 
rsm0  r  •  r 

where  <7  =  |  g  |  and  p  is  the  density. 

For  the  plasticity  model,  the  von  Mises  yield  condition  is  assumed  to  hold.  Ex¬ 
pressed  in  terms  of  the  principal  stresses,  ait  i=l,2,3,  i.e.,  the  eigenvalues  of  T,  this 
condition  reads 

(2.2)  (<n  -  p)2  +  («r2  -  p)2  +  (its  ~  P)2  =  2 p2s2, 
where 

P  =  ^trT  =  ^(^rr  +  +?W)  =  -(Cl  +  <72  +  ^3) 

is  the  average  stress  and  s  =  sin  5,  <5  being  the  angle  of  internal  friction. 

A  flow  rule  completes  the  model.  The  eigenvectors  of  the  strain  rate  tensor  V 
and  the  deviatoric  part  of  the  stress  tensor  T  are  assumed  to  be  parallel.  This  is 
the  Levy  flow  rule,  which  can  be  equivalently  expressed  as  the  existence  of  a  positive 
scalar  function  A  >  0  such  that 

(2.3)  V  =  X(T-pI). 

The  alignment  condition  of  the  eigenvectors  of  T  and  V  in  effect  neglects  the 
rotation  of  a  material  element  during  deformation,  a  controversial  assumption.  There 
is  experimental  evidence  that  misalignment  may  occur  under  some  circumstances. 
Alternative  models  which  allow  for  the  above  eigenvectors  to  be  somewhat  out  of 
alignment  have  been  proposed,  see  e.g.  [20].  However,  to  the  best  of  our  knowledge, 
there  does  not  seem  to  be  enough  experimental  data  to  favor  one  type  of  models  over 
the  other.  We  refer  to  [7]  for  a  lucid,  if  somewhat  dated,  account  of  the  situation. 
Further,  and  more  importantly,  the  fact  that  we  look  here  exclusively  at  radial  flows, 
renders  the  distinction  between  the  two  types  of  model  less  of  an  issue.  Indeed,  simple 
explicit  calculations  show  that  in  the  case  of  radial  flows  in  three  dimensional  conical 
hoppers  of  circular  cross  section,  for  instance,  the  misalignment  predicted  by  Spencer’s 
double  shearing  model  [20]  is  less  than  10°  in  all  cases  of  physical  interest. 

The  unknowns  characterizing  the  flow  are  the  velocity  v  and  the  stress  tensor 
T,  assuming  constant  density.  They  are  determined  by  (2.1)  which  corresponds  to 
conservation  of  momentum,  the  yield  condition  (2.2)  and  the  above  flow  rule  (2.3). 
Note  that  (2.3)  also  implies  incompressibility,  and  thus  conservation  of  mass,  as 

V  •  v  =  —try  =  Atr(T  —  pi)  =  0. 


The  type  of  the  system  (2. 1-2. 3),  which  has  to  be  supplemented  with  side  conditions 
discussed  below,  can  be  determined.  In  [17],  it  was  shown  that,  depending  among 
other  things  on  the  internal  friction,  the  system  for  fully  three-dimensional  flows  can 
be  mixed  hyperbolic-elliptic,  elliptic  or  have  no  definite  type. 


4 


P.A.  GREMAUD 


Following  Jenike  [8],  [9],  we  seek  similarity  solutions  with  radial  symmetry .  More 
precisely,  the  velocity  field  is  assumed  to  be  purely  radial,  i.e. 

[Vr.Wfli^]  =  far, 0,0], 

while  the  stress  components  are  of  the  form 


Trr  —  rTrr($,0), 

and  similarly  for  the  other  stress  components.  The  strain  rate  tensor  simplifies  to 


V  = 


drVT  2r^Vr  2rs\n6®'PVr 


“£®r 


0 


-rtV 


The  flow  rule  (2.3)  has  several  fundamental  consequences.  First,  since  again  V-u  =  0, 
we  obtain 

drvr  +  -vr  =  0. 


The  proper  scaling  for  the  velocity  is  consequently 

vr{r,6,<t>)  =  -\v(9,<f>), 


where  we  further  have 
(2.4) 

Finally,  (2.3)  implies 


dgV  = 


2  Tre 

— - v. 

Tgg  -  p 


Tgg  —  *^"00  and  Tg^  0. 


The  yield  condition  (2.2)  can  be  expressed  in  terms  of  the  four  stress  unknowns  rrr, 
Tr0 ,  rrip  and  Too  by  using  the  invariants  of  T,  see  (2.8)  below.  This  results  in  the 
following  system  of  three  partial  differential  equations  and  one  algebraic  constraint  in 
the  four  remaining  stress  unknowns  rrr,  Trg,  rr<j  and  Tgg 


(2.5) 

(2.6) 

(2.7) 

(2.8) 


— +  dgTrg  +  3rrr  -  2 Tgg  +  Trg  COt  6  =  -pg  COS  6, 
sin  9 

-r-rd^Tgg  +  4rr^  =  0, 

sin  9 

dgTgg  +  4  rrg  =  pgsln9, 


(Trr  —  Tgg)2  +  3  Tr0  +  3  T, 


=  ifar 


+  2  Tqq)^  52 


The  yield  condition  (2.8) 
TtAT=-Zt^ 


deserves  further  comments.  It  can  be  written 


Trr 

r  1-^ 

1  3 

0 

_1  _  2*1  I 

1  3 

TVQ 

,  a  = 

0 

3 

0 

T 80 

1  2s2 

L  1  3 

0 

1  4  s2 

1  3  J 

with  T  = 
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Fig.  2.1.  Two  representative  level  surfaces  corresponding  to  the  yield  condition.  In  each  case , 
only  the  physically  relevant  nappe  is  showed. 


The  symmetric  matrix  A  can  be  diagonalized 


'  Ax 

‘ 

A  —  XAXt,  with  A  = 

a2 

A3  . 

,  *  = 

Xx 

X2 

*3 

where 

A1=3,  A2,3  =  l-^-±^/l  +  |s2  +  (^-)  , 

—  [0, 1,0]T,  *2,3  =  (jY  ±  |>/36  +  48s2  +  25s*)  ,0,  l]  . 

As  is  easily  seen  from  those  expressions,  Ai  and  A2  are  positive  for  any  S  >  0,  while 
A3  is  negative.  Let  S  =  [  £  77  £  ]T  =  XTT\  the  yield  condition  takes  the  form 

(2.9)  StA  S  =  3£2  +  A2t f  +  A3C2  =  -3  r2^, 

which  corresponds  to  a  family  of  hyperboloids  of  two  sheets  parametrized  by  |rr^|, 
see  Figure  2.1.  The  case  Tr<p  =  0  is  a  cone.  Only  one  of  the  two  sheets  is  physically 
relevant.  Indeed,  granular  materials  can  only  support  compressive  stresses,  i.e.,  Oi  > 
0,  i  =  1,2,3.  This  condition  is  satisfied  on  only  one  of  the  sheets,  provided  S  <  60°. 
The  corresponding  result  given  below  can  be  found  in  [17]  where  it  is  stated  and 
proved  using  exclusively  the  properties  of  the  principal  stresses  (i.e.,  (2.2)  instead  of 
(2.8)). 

Lemma  2.1.  For  any  value  of  r^,  the  hyperboloid  of  two  sheets  (2.8)  has  one 
sheet  corresponding  to  compressive  stresses ,  Oi  >  0,  i  =  1,2,3,  if  and  only  if  the 
internal  angle  of  friction  satisfies  0°  <  S  <  60°. 
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Proof.  As  a  direct  consequence  of  the  structure  of  the  scaled  stress  tensor  T  = 

Trr  TtQ  Tr^> 

rre  Tee  0  ,  it  can  be  seen  that  the  principal  stresses  are 

Tr<f>  0  Tqq 

(2.10)  d\  =700, 

(2.11)  ^2,3  =  —  y  —  ±  \J (t -09  -  rrr)2  +  4rr2fl  +  4rr20, 

without  assuming  any  ordering  of  the  <5Vs.  The  component  has  thus  to  be  positive. 
Further,  we  have 

detT  —  tqq[ttttqq  —  r2#  —  t2^)  —  5*i5,2^3i 

and  thus 

(2.12)  Trrre0  -  =  0-2CT3, 

Therefore,  cr2  and  03  only  have  the  same  sign  if  rrr  >  0.  Solving  for  rrr  directly 
from  the  yield  condition  (2.8)  shows  that  rrr  >  0  for  positive  values  of  tqq  provided 
2s2  -  3\/3 s  +  3  >  0,  see  (2.13)  below.  Since  s  —  sin  5,  it  is  equivalent  to  S  <  60°. 

Combining  (2.8)  with  (2.11),  we  observe  that  the  remaining  b<i  and  <73  are  positive 
if 

c2  4s2  452 

(y  -  1  )Trr  +  (y  -  l)Tge  +  (-3-  +  2 )ttttbo  <  3rrrTfi9, 

which  is  clearly  satisfied  if  again  8  <  60°.  □ 

At  this  point,  the  system  (2. 5-2.8)  can  be  rewritten  in  several  ways.  One  could  for 
instance  introduce  a  la  Sokolovskii  variables  [14]  which  would  essentially  correspond 
to  solving  (2.9)  by  replacing  £  and  77  by  one  new  variable 

f  =  -yyC^C2  +  3r^)  sin  2-0,  V  =  yy^sC^  +  Sr^)  cos  2ip. 

Note  that  the  above  expressions  make  sense  because  of  the  yield  condition  (2.9)  and 
the  fact  that  A2  >  0.  The  approach  would  lead  to  an  elliptic  system  of  three  fully 
nonlinear  first  order  equations  in  the  variables  £,  tp  and  rr^.  We  choose  rather  to 
work  with  the  original  scaled  variables  after  having  solved  (2.8)  for  rrr 

7 ~rr  =  {TtO  j 

(2.13)  =  — yj  ((3  +  2s2)t99  ±  3^3s2r2fl  -  (3  -  s2)rr2fl  -  (3  -  s2)^)  . 

As  can  obviously  be  seen  from  both  the  above  expression  and  Figure  2.1,  for  any  triple 
(Tr0,Tr<£,T00),  one  can  have  zero,  one  or  two  corresponding  values  of  rrr.  The  zero 
solution  case,  i.e. 

rfe  <  yr(Trf  + 

corresponds,  by  assumption  here,  to  stress  states  incompatible  with  the  physical  prop¬ 
erties  of  the  material  under  consideration. 
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The  two  solutions  in  (2.13)  are  related  to  the  so-called  active  (+  sign)  and  passive 
(-  sign)  states  of  the  flowing  material,  see  e.g.  [14].  As  the  passive  state  is  the  one 
that  tends  to  be  observed  experimentally  upon  discharge  of  the  hopper,  we  adopt,  to 
fix  the  ideas,  rTT  =  /_(rr0,Tr0,T00).  However,  the  analysis  presented  below  applies 
equally  to  both  cases. 

The  unknowns  rT$  and  ttq  can  be  eliminated  from  the  system  by  using  (2.6)  and 
(2.7).  Relation  (2.5)  becomes 

.  "T'A^4><t>TQ0 - : — -deismOdeToe)  H-  12rrr  —  8700  =  —6pg  cos  6. 

sm  6  sin  v 

One  recognizes  the  above  differential  operator  -^4^ (sin Ode-)  =  -A  as 
Laplace’s  operator  on  the  sphere.  The  system  is  closed  by  expressing  rTT  as  a  function 
of  700  and  its  first  derivatives  through  (2.6),  (2.7)  and  (2.13),  leading  to 

(2.14)  —  Ar00  +  xree  =  -6pg  cos  <9  +  F(0,  ree^doree,  d^reo), 

where  x  ~  4-^f/a-  and 


F(0,T00,50T00,90T00)  -  3^2  y 48 s*r*e  -  (3  -  S2)  ((pg sine  -  deree)2  +  (^)2)  • 

The  above  nonlinear  Helmholtz  equation  (2.14)  has  to  be  solved  in  the  domain 
ft fuu  corresponding  to  the  intersection  of  the  unit  sphere  centered  at  the  vertex  of 
the  hopper  with  the  hopper  itself.  Whenever  possible,  symmetry  properties  are  taken 
into  account.  The  resulting  computational  domain  ft  is  a  “spherical  triangle”,  with 
two  edges  corresponding  to  arcs  of  great  circles  and  a  third  one  corresponding  to  the 
outer  boundary,  see  Figure  2.2  in  the  case  of  a  domain  with  rectangular  cross  section. 
The  boundary  9ft  of  ft  consists  of  the  four  parts 

9ft  ~  To  u  Ti  ur2  ur3, 

where  T0  is  the  trivial  boundary  {(0,0); 6  =  0}  and  rij3  =  {(0,0); 0  <  9  <  9W,4>  = 
T0u;},  with  6W  >  0  and  <j)w  >  0.  The  boundary  T2  is  represented  as  the  graph  of  an 
even  function  C  of  class  C1,  9W  =  C(±<pw) 

r2  =  {(0,<j>)-,6  =  C(4>),~<pw  <4><  4>w). 

The  domain  ft  in  which  the  problem  is  to  be  solved  is  thus 

ft  =  {(M);  ~<t>w  <  <t>  <  K,0  <  6  <  C{<t>)}, 


see  Figure  2.2. 

The  boundary  conditions  are  derived  from  physical  considerations.  On  Ti  and 
r3,  we  impose 

(2.15)  T00(*,  0-u/)  T## (' j  0tt;) 5  d<f}'T0Q(<')  0it/) —  O(pT00  (*,  0-u;)  “  0. 

Indeed,  by  symmetry,  it  is  clear  that  one  should  have 


T0d{’)  4*w)  — ^0^(*7  0ty))  ^0T00(*,  (f)w )  —  0u;)* 
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Fig.  2.2.  Geometry  of  the  computational  domain. 


Further,  by  a  formal  regularity  argument,  one  obtains  that  d^ree  vanishes  on  Ti  and 
r3.  Indeed,  if  one  solves  the  problem  in  ftfuii  instead  of  ft  and  freezes  the  coefficients 
in  the  nonlinear  part  F  of  (2.14),  classical  regularity  results  apply,  see  [6]  §3.2  and 
Remark  3.2.4.6.  This  yields  iJ^-regularity  of  the  solutions.  For  such  a  solution  it,  one 
has  on  each  4>  =  cst  line  T# 

=  700iT, 

where  correspond  to  the  values  from  each  side  of  T#  and  the  mapping  7  :  H2 (ft juu)  ~> 
Jff1/2(r^)  which  takes  u  to  7 d$u  is  the  normal  trace  operator,  see  [6],  Lemma  1.5. 1.8. 

The  condition  d(j>Tee(±^>w)  =  0  has  also  been  verified  numerically  by  varying  the 
domain  of  resolution  (for  instance,  by  rotating  it  by  7r/4  in  case  of  a  square  domain). 
On  T21  the  law  of  sliding  friction  applies,  i.e. 

|  Tt  I  =  PwTjy,  on  r2, 

where  f  t  and  Tn  are  the  scaled  tangential  and  normal  stresses  on  the  hopper  wall, 
and  fiw  is  the  coefficient  of  wall  friction.  Since  the  outer  unit  normal  to  the  wall  is 

N  _  [0,  sin  0, —£'(<£)] 

y^sin  20  +  C'(0)2 

it  follows  that 

Tn  =  Tee  and  Tr  =  WeTre  4-  N^t^,  0,0]T. 

The  boundary  condition  on  r2  then  reads  by  (2.6),  (2.7) 

pg  sin2  C(<j>)  -  sinC(4>)deT0e  +  d^Tee  = 

(2.16)  -4 pw  \/sin2C(</>)  +C'(0)2tm, 

which  is  essentially  a  Robin  boundary  condition. 
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3.  Mathematical  analysis.  The  above  problem  presents  several  nonstandard 
features  and  difficulties  (Laplace-Beltrami  operator,  Robin  boundary  condition,  non 
Lipschitz  coefficients,  restricted  range  of  admissible  values).  In  this  section,  an  exis¬ 
tence  and  local  uniqueness  result  is  established. 

We  equip  L2(Cl)  with  the  inner  product  (u,v)  =  j£W  uv  sin  9  ddd(p.  Let 

V  =  {v  6  Vu  =  [devt  €  L2{ fi)2, 

u  (•>  ~<t>w)  =  u(-,  <pw)  in  the  sense  of  the  trace}. 

The  space  V  is  equipped  with  inner  product  and  norm 


(u,v)v  =.(u,v)  +  (Vu,Vv),  IMIv  =  v'(u,u)v. 

Formally,  if  u  stands  for  a  classical  solution  to  problem  (2.14)  and  if  v  €  V,  we 
obtain  after  integration  by  parts  and  use  of  the  boundary  conditions  (2.15)  and  (2.16) 

rr  r<f>w  rC(<f>) 

-J J  Au  v  dui  =  J  J  Vu  •  Vv  sin  6  d6d<f> 

~  /  0  (4/X“  \/sin2C^)  +  C'(0)2u(C(0),  <f)  +  PPsin2  Ci^viCicf)),  <j>)  d<j>. 
Let  a  :  V  x  V  -»  R  be  the  bilinear  form  defined  by 

/4>w  pC((f)) 

/  {  Vu  •  Vu  +  xuv)  sin  0  d9d(p 

-0U >  d o 

/<pw  j - 

V sin2  C(4>)  +  C'{4>Y  u(C(<t>),  <t>)v(C{4>),<t>)  #. 

For  T2,  i.e.,  C  of  class  C1,  it  is  well  known  that  the  trace  operator  7  :  V  iJ1/2(r2) 
is  linear  and  continuous  [6],  Th.  1.5. 1.3,  p.38.  Consequently,  the  form  a  :  V  x  V  R 
is  continuous.  Further,  under  a  condition  that  links  wall  friction,  internal  friction  and 
the  geometry  of  the  domain  (a  precise  condition  is  given  as  part  of  the  proof  of  the 
next  result),  the  bilinear  form  a(-,  •)  is  V'-elliptic. 

Proposition  3.1.  The  bilinear  form  a(v)  is  V-elliptic,  provided  T2  is  a  curve 
of  class  C 1  and  is  sufficiently  small. 

Proof  Let  e$  =  [1,  0]  and  h{(j>)  —  sin2  C(0)  -1-  C'(0)2;  we  have 

JJ h(<fi)  V(u2)  *  eo-du  =  2  J j h{(j))udeud(jj. 
n  q 

On  the  other  hand,  the  divergence  theorem  leads  to 

ff  h(4>)  V(u2)  •  eg  dui  =  -  j  j  u2  V  •  ee)du  +  J  h(^)u2ee  ■  N  da 


=  J  h(<j>) 

r2 


u2ee  •  N  da , 
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where  N  stands  again  for  the  unit  outer  normal  and  where  we  have  used  (2.15).  On 
T2,  we  observe  that 


h{cj>)  eg  •  N  =  ^/sin2  C(</>)  +  C'(<£)2  eg  ■  N  =  sin C(<f>),  <f>  €  [~<j>w,  </>w], 

and  thus 

f  h(<j>)u2eg  ■  N  dcr  >  ■  - — _sm  ®nun  f  h{<j>)u2da, 

r2  y  S*n  ®ma.x  +  £*00  r2 

where  6min{max)  =  min^€[_0u;>^](max)C(^)  and  C ^  =  max0€[_0wt0iii]  |  C'(0)  |.  Com¬ 
bining  the  last  three  expressions  and  using  da  =  sin  6  d(p  yields,  for  any  rj  >  0 


\/sin2  ^mai  +  c^2  rr 

/  h(<l>)u2(C(<t>),<t>)d<P<  2+ - -f- - 

J  sin  emin  JJ 


|h(0)|  |u|  dut 


<  Sill2  dmax  +  C'^  (  1 


sin2  9n 


^Nli*(0)  +  v\\dgu\\lHQ)y 


From  the  definition  of  a(-,  •),  sufficient  conditions  of  F-ellipticity  are  then  found  to  be 


X- 


4/iy;  sin  Ornax  *t* 


>0  1  -  i{2wri 


sin  ^TTiax  C-r 


V  smzemin  "  "  'r'u/w  sin2 


>0, 


for  some  rj  >  0.  A  sufficient  condition  is  then  for  instance 

,2  a  .  nt  2  \  2 


X-32A(^?;+C”  )  >0. 
\  sin  Ojyiin  J 


We  can  now  state  a  weak  formulation  of  our  problem:  Find  a  function  u  €  V  such 

that 

a(u,t;)  =  JJ F(u)vdu  -6pgJJ  cos6vdu  +  pg  J  sin2  C(</>)v(C (<{>),  (f>)  d(j>, 
n  Q 

for  any  and  where,  with  a  slight  abuse  of  notation,  F  still  denotes  the  nonlinear 

function  (2.14).  We  define  A  £  £(K,  V’*)  and  $  £  K*  by  respectively 

<  Au,v  >  =  a(u,v),  Vu,i;  £  V, 

<*'”>— 6  PSJI  cos  da?  +  /:  sin2  C(0)t?(C(0),  0)  d(j> ,  Vu  £  V, 

n  ^ 

where  <  *,  •  >  denotes  the  V",  F*-duality  product.  Identifying  in  the  obvious  way  F(w) 
with  an  element  of  V*,  the  problem  becomes:  Find  u  £  V  such  that 

(3.1)  Au  =  F(u)  +  $. 

Let  us  denote  b  y  Q  :V  ^  V*  the  mapping 


G(u)  —  Au  —  F(u) 
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* 

and  let  u0  =  -pg  cos  8  e  V.  Note  that  Q{u0)  €  V*  is  well  defined  as  F(u0)  is  itself 
well  defined  by  construction,  provided  0max  <  7r/2.  Direct  calculations  show 

ay,v*)  BDg(u)  =  a-df(u), 


where 


<  DF(u)h,v  >= 


9  [  f  48s2  uh- (3 -s2)[-(pg sm6 -deu)deh  +  &g£]  , 

2  _  ^2  II  I  -  ^ 

n  y48 s2u2  -  (3  -  s2)[(pgsmO  -  deu )2  + 


for  any  h,v  £  V. 

PROPOSITION  3.2.  Under  the  assumptions  of  Proposition  3.1,  there  are  neighbor¬ 
hoods  U  C  V  and  U *  C  V *  of  uo  and  G{uq)  respectively  such  that  for  each  $  €  U*, 
the  functional  equation  (3.1)  has  one  and  only  one  solution  in  U. 

Proof.  Prom  the  previous  expression  and  the  Inverse  Mapping  Theorem,  Q  is  a 
C^diffeomorphism  between  a  neighborhood  of  u0  in  V  and  one  of  G(uo)  in  V*.  The 
proposition  follows.  □ 

Arguments  following  the  lines  of  the  formal  remarks  at  the  end  of  §3  could  be 
used  to  obtain  #2(fi)-regularity  of  the  solutions.  We  do  not  pursue  this  issue  further 
here. 

4.  Numerical  analysis.  In  order  to  simplify  the  numerics,  the  problem  is 
mapped  onto  a  rectangular  computational  domain.  Prom  now  on,  the  function  C 
which  describes  T2  is  assumed  to  be  of  class  C2.  We  define  the  new  coordinates 


e  =  ew 


e 

m 


and  <£  =  (/>. 


Note  that  {r,  0,$}  is  not  an  orthogonal  coordinate  system.  Keeping  in  mind  that 
6  =  QC{4>)/dw1  the  problem  for  the  transformed  unknown  £/(©,$)  =  u{Q,<j))  takes 
the  form 


~d**U  +  2Q^de*U  -  (02C'($)2  +  02  sin2  6)  deQU 

+  (c^)  [c"($)c($)  “ 2 C'($)21  ~  \  d~Y^y-)  9®u  +  x sifl2  0  U 

(4.1)  =  -6pg cos 9 sin2  9  +  F(Q,  <f>,U,d@U,  d<s>U),  (0,  $)  e  (0,6W)  x  <f>w), 


where 

9  sin  6 
*  ~  3-s2 

The  boundary  condition  on  1?2  becomes 
pgsin2C{§) -sinC($)  f^dQU  +  C  ® 


48s2  sin2  9  U2  -  (3  -  s2)  (pg  sin2  9  -  d@U]2  +  [d*U  -  ^~-d@U]2^ . 


cm 


sinC($) 


a  cm, 

“w  oq  +  0$ 


(4.2) 


-4  sm2  cm +  C'm2u,  §  e  {-KAw). 


The  above  problem  is  discretized  by  collocation;  Chebyshev  collocation  at  the 
Chebyshev-Gauss-Radau  points  in  used  in  0,  while  Fourier-cosine  collocation  at  the 
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Fourier  collocation  points  is  used  in  <I>.  More  precisely,  we  set 

N- 1  M/2-1 

(4.3)  CW®,*)=E  E  UnmMQ)eim{i*+nK 

n= 0  m=  —  M f2 


where  {^n}n= o  are  Lagrange  interpolation  polynomials  at  the  Chebyshev-Gauss- 
Radau  nodes  on  [0,0^],  i.e. 


(4.4) 


H-  cos( 


27 xj 
2N  — 


j  =  0, . . .  ,iV  -  1. 


This  choice,  as  opposed  to  the  more  standard  Chebyshev-Gauss-Lobatto  collocation, 
see  e.g.  [2],  §2.4,  results  from  the  nature  of  the  boundary  condition  along  0  =  6W. 
For  completeness,  we  derive  the  expression  of  the  collocation  derivative  below  (which 
we  have  not  been  able  to  find  in  the  literature). 

Lemma  4.1.  The  Lagrange  interpolation  polynomials  on  the  Chebyshev- Gauss- 
Radau  nodes  (44)  ^re  9iven  h 


(ff *  (f  -  0  +  ITT Tr»-  (£  -  0) '  ’  =  ° N-h 


where 


co  =  1  -  2  AT, 


ew  /  2tt Nj  ...  ..  2ir(N  —  l)j 

CJ  =  ~2§-  [Ncos  +  (N-l)  COS 


2JV-1 


2  N  -  1 


where  T^(x)  =  cos(iVarccosx),  |x|  <  1,  is  the  Chebyshev  polynomial  of  degree  N. 

The  above  result  can  easily  be  verified  through  the  use  of  l’Hospital’s  rule  and 
elementary  properties  of  the  Chebyshev  polynomials.  Interpolation  at  the  nodes  (4.4) 
of  a  function  u  of  0  defined  in  [(),#,„]  simplify  takes  the  form 


JV-l 

iNu(&)  =  Eu(0^(0)- 

j= 0 

By  definition,  the  Chebyshev  collocation  derivative  of  u  at  those  nodes  is  then 


N- 1  N- 1 

(W(0/)  =  E  u(QiWQ‘)  =  E  VIA®/), 

j= 0  j=0 

with  Vij  =  x/;'(0/).  The  collocation  derivative  at  the  nodes  can  then  be  obtained 
through  matrix  multiplication.  Elementary  albeit  tedious  calculations  lead  to  the 
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following  expressions 

f  hh^N-i), 

1  _ 2 


vij  =  { 


iU  =  j  =  0, 


cq&u,  sin2  “i 


Cj{_  1 

Cj  &  i  Q  j 


IT  cos  MiL  +  (N-1)  cos  if  j  =  0, 

1  =  1, . . . ,  TV  —  1, 

ifj  =  l,...,TV-l, 
j  5^  Z  =  0, . . .  j  N  —  1, 


1  3flu,  — 2Qj  flu, l 

4 ©j  —  ©j  4  c  " 


(iV  —  l)2  sin((Ar  -  1)  arccos(^-  -  1))  j 


(^sinCTVarccos -  1))+ 


if  /  =  T  =  1, . . .  ,TV  -  1. 


In  the  $  direction,  the  collocation  points  are  taken  as  the  usual  Fourier  collocation 
nodes,  i.e., 


21 

(4.5)  l  —  0, . . . ,  M  -  1. 

Let  U  be  the  TV  x  M  matrix  of  coefficients  Unm,  n  =  0, . . . ,  TV  -  1,  m  =  0, . . . ,  M/2  - 
1,  —M/2, . . . ,  —  1  and  let  W  be  the  M  x  M  Fourier  matrix 


"  1 

1 

1 

..  i 

1 

Um 

•  ■  “S-1 

1 

UM 

-r-' 

_  1 

2(M— 1) 

Urri  '  . 

where  u>m  =  e'2jr/M  is  the  primitive  M-th  root  of  unity.  Further,  if  C  is  the  M  x  M 
diagonal  matrix  with  diagonal  [0, . . . ,  M/2  -  1,  -M/2, . . . ,  -1],  then  for  any  j,l,  j  = 
0, . . . ,  TV  —  1,  l  =  0, . . . ,  M  —  1,  the  nodal  values  of  Unm  and  its  derivatives  can  be 
expressed  as  follows 


UNM{®j,$l)  =  (UW)jl, 
d$UNM(®j,  $()  =  4i  (UCW)ji, 
d<s>$UNM(Qj,$i)  =  —16  (UC2W)ji, 
deUNM{®j^i)  =  (VUW)jh 
de$UNM{®j,  $i)  =  4i  ( VU£W)ji , 
deeUNMiQj^t)  =  (V'2UW)ji. 

The  TV  x  M  matrix  of  unknown  coefficients  Tf  clearly  satisfies  a  matrix  equation 
of  the  type 


AiUBi  +  . . .  ApUBp  =  F{U), 

where  the  A/s  are  TV  x  TV  matrices  while  the  T?;’s  are  M  x  M.  The  above  system 
is  obtained  by  enforcing  the  conditions  that  first,  the  discrete  solution  UNM  from 
(4.3)  satisfies  the  PDE  (4.1)  at  the  collocation  points  {(Qj, $/)},  j  =  1,...,TV  -  1, 
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l  =  0, . . . ,  M—  1  and  second,  that  the  boundary  condition  (4.2)  is  verified  at  the  nodes 
{0o,  $/)},  /  =  0, . . . ,  M  —  1.  Note  that  no  side  conditions  of  any  form  are  imposed  in 
the  neighborhood  at  ©  =  0.  We  denote  by  A<g)B  the  NMxNM  matrix  corresponding 
to  the  Kronecker  product  of  an  TV  x  TV  matrix  A  by  an  M  x  M  matrix  B .  Further,  for 
any  matrix  A ,  vec(A)  denotes  the  vector  formed  by  stacking  the  columns  of  A.  One 
can  check  [12],  p.410 

vec (AUB)  =  ( Bt  ®  A)  vec (W). 

A  direct  consequence  of  this  elementary  relation  is  that  the  above  matrix  equation 
can  be  rewritten 

(4.6)  HI  vec  (U)  =  vec (F(W)), 

where  H  =  ®  A/)  is  a  non  sparse  7VM  x  iVM-matrix. 

The  above  nonlinear  system  is  numerically  solved  as  follows.  First,  the  matrix 
HI  is  factorized  into  HI  =  LV,  where  V  is  an  upper  triangular  matrix  and  L  is  a 
“psychologically”  lower  triangular  matrix  (LU  factorization).  Then,  the  nonlinear 
equation 

vec(li)  =V~1L~1vec(F(U)), 

is  solved  by  a  Newton-GMRES  solver  [10]. 

5.  Numerical  results.  The  numerical  approach  is  tested  by  comparing  with  an 
exact  solution  obtained  for  a  simplified  linear  problem  in  a  domain  corresponding  to 
a  circular  cone.  More  precisely,  we  consider  the  following  problem 

(5.1)  _  sin^(sin  6u,(6)Y  +  xu  =  F  for  0  <6  <8W, 

(5.2)  u7(0)  =  0  u'{0w)  -  pg sin 6W  +4 pwu(6w), 

where  F  is  taken  as  constant.  The  above  equation  and  boundary  conditions  of  this 
model  problem  corresponds  exactly  to  (4.1,  4.2)  whith  a  very  simplified  right-hand 
side. 

The  change  of  variable  v(z)  =  u(8)  with  z  =  cos0  leads  to 
((1  -  z2)v'(z))’  +  v{v  +  1  )v(z)  =  -F, 

where  v(v+ 1)  =  — x*  A  fundamental  system  of  solution  to  the  homogeneous  equation 

((1  -  z2)v'(z)y  +  v[y  +  1  )v{z)  =  0, 

is  provided  by  the  Legendre  functions  Pt/(z)  and  Qu(z)  of  the  first  and  second  kind 
respectively,  see  e.g.  [3],  §8.82,  8.83.  The  coefficient  v  defined  by  v(v  +  1)  =  ~x 
satisfies 


where  A  —  yjx  ”  \  a  real  parameter  as,  obviously,  x  =  4 >  4  for  any  material. 
We  can  thus  get  the  solution  in  terms  of  P_i+iX(z)  and  Q_x+iX(z)  which  are  conical 
functions ,  see  [3],  §8.84.  The  general  solution  to  (5.1)  is  then 


(5.3) 


u{9)  -  AP_i_+iX( cosd)  4-  BQ_i+iX(cos9)  +  F/X 


FLOWING  GRANULAR  MATERIALS 


15 


N  =  M 

8 

10 

12 

14 

error 

7.061  (-8) 

1.723  (-10) 

3.567  (-13) 

5.995  (-15) 

Table  5.1 

Maximum  norm  of  the  error  for  the  model  problem  (5.1),  (5.2). 


where  the  last  term  is  a  particular  solution  to  the  considered  problem.  Since  for  real 
6 ,  Q_i+iA(0)  is  always  imaginary  while  P_i+iA  is  always  real,  we  get  B  =  0.  After 
noticing  that  u'(0)  =  0,  the  value  of  A  can  be  found  from  (5.2) 


A  _  pgsin6w  +  4:{j,wF/x 

pl  tw) 


where  p(6)  =  P_i+iA(cos0). 


Our  numerical  solver  is  now  applied  to  the  present  model  problem  (5.1),  (5.2). 
The  corresponding  numerical  solutions  are  compared  to  the  exact  solution  (5.3).  The 
results  are  summarized  in  Table  5.1,  where  the  maximum  norm  of  the  error  is  given 
as  a  function  of  the  number  of  nodes,  taken  here  as  N  =  M.  This  illustrates  two 
facts.  First,  as  expected,  the  method  is  found  to  be  spectrally  accurate.  Second,  for 
this  simple  problem,  a  mesh  as  small  as  14  x  14  leads  to  errors  of  the  order  of  the 
round-off  errors. 

The  standard  fully  axisymmetric  case  of  a  circular  cone  for  the  full  problem 
provides  an  additional  way  of  testing  the  approach,  this  time  in  the  nonlinear  regime. 
Again,  axisymmetry  considerably  simplifies  the  problem.  Since  no  variable  depends 
on  (j>  in  that  case,  a  direct  consequence  of  the  flow  rule  (2.3)  is  that  the  (scaled)  stress 
tensor  T  reduces  to 


f  = 


T~rr  l~r8  0 

Trd  Tqq  0 

0  0  Tqq 


Further,  the  yield  condition  (2.8)  reduces,  in  (rrr, ttq,tqq )-space,  to  the  cone  illus¬ 
trated  in  Figure  2.1  ( T ^  =  0),  as  opposed  to  the  family  of  hyperboloids  discussed  in 
§2.  This  cone  can  be  parametrized  by  the  Sokolovskii  variables  (cr,V ;)  where  a  is  the 
average  stress  and  ^  is  a  new  variable,  as  follows  [9],  [14],  [17] 


2  1 
rrr  =  cr(l - -j=s  cos  2^),  rr0  =  -as  sin  2^,  tqq  =  a(l  +  -=s  cos  2*0). 


vr 


The  above  parametrization  “solves”  the  yield  condition.  The  remaining  unknowns 
(a,^)  are  determined  by  solving  the  equations  of  equilibrium  (2.1),  which  now  only 
yield  two  ordinary  differential  equations.  The  corresponding  boundary  value  problem 
can  easily  be  solved  numerically.  We  omit  the  details,  see,  e.g.,  [5],  A  comparison  in 
terms  of  the  average  stress  a  between  results  from  the  ODE  code  using  the  Sokolovskii 
variables  on  the  one  hand,  and  results  from  the  present  approach  on  the  other  hand, 
is  illustrated  in  Figure  5.1.  The  agreement  is  excellent  (relative  error  <  .02%). 

Having  gained  some  confidence  in  the  approach,  we  now  apply  it  to  the  general 
full  problem.  Results  are  presented  for  the  following  kind  of  hoppers.  For  0  <  A  <  1, 
we  consider  the  following  family  of  domains 

C(4>)  =  (1  -  A )0W  +  A  arctan( tan  —  C°S  — ) , 

COS  (p 
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Fig.  5.1.  Comparison  between  ODE  solver  (solid  curve)  and  present  spectral  PDE  solver  (o) 
for  a  circular  (axisymmetric)  cone ;  9W  =  30°,  S  =  30°,  pw  =  tan  (15°),  N  =  M  =  20. 


corresponding  to 

•  conical  hopper  with  circular  cross-section  A  =  0, 

•  pyramidal  hopper  with  transitional  cross-section  0  <  A  <  1, 

•  pyramidal  hopper  with  square  cross-section  A  =  1, 

In  each  case  treated  below,  <f>w  =  7r/4  and  6W  is  a  free  parameter,  both  angles  being 
defined  in  Figure  2.2. 

The  results  presented  below  are  meant  to  illustrate  the  qualitative  effects  of  the 
geometry  of  the  hopper  on  the  flows.  All  results  are  represented  on  horizontal  cross- 
sections  taken  at  the  same  vertical  height.  This  requires  appropriate  scaling  according 
to  the  principle  outlined  in  §2  (stress  =  O(r),  velocity  =  0(r~2)).  Unlike  the  stress, 
the  velocity  field,  which  is  computed  from  (2.4),  is  only  defined  up  to  a  scalar  multiple. 
This  indeterminacy  can  be  solved  by  imposing  an  outflow  rate  for  instance.  Inciden¬ 
tally,  in  practice,  the  flow  through  outlets  at  the  bottom  of  the  hopper  is  almost 
always  imposed,  through  some  kind  of  feeder  device  for  instance.  Here,  the  velocity 
fields  are  normalized  to  have  maximum  value  1.  In  Figure  5.2,  a  comparison  between 
a  circular  conical  hopper,  a  transitional  hopper  A  =  1/2  and  a  square  pyramidal  one 
with  “same”  opening  angle  is  offered.  By  same  angle,  we  mean  that  at  a  given  com¬ 
mon  height,  the  geometric  domains  corresponding  to  horizontal  cross-sections  admit 
inscribed  circles  of  equal  diameter.  The  material  constants  are  chosen  as  6  =  40°  and 
fiw  =  tan(20°). 

All  calculations  were  done  on  a  quarter  domain  (the  solid  diagonal  line  appearing 
in  all  figures  is  an  artifact  of  the  plotting  software).  Several  observations  can  be  drawn 
from  the  above  experiment.  First,  flows  in  pyramidal  hoppers  seem  to  exhibit  larger 
stresses  than  flows  in  purely  conical  containers.  To  derive  more  precise  practical  con¬ 
clusions,  arguments  related  to  the  volume  of  material  effectively  contained  would  have 
to  be  considered.  Second,  and  somewhat  surprisingly,  the  presence  of  corners  seems  to 
have  relatively  little  effect  on  the  velocity  flow.  Comparisons  with  both  laboratory  ex¬ 
periments  and  different  plasticity  models  would  be  extremely  interesting.  In  the  usual 
axisymmetric  case,  the  latter  kind  of  comparison  can  be  easily  performed.  Models 
based  on  a  TYesca  yield  condition  for  instance,  rather  than  the  present  von  Mises  con- 


average  stress 


velocity 


-0.5  -0.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3  0.4  OS 


-0.5  -0.4  -0.3  -0.2  -0.1  0  0.1  0  2  0.3  0.4  O.S 


Fig.  5.2.  Horizontal  cross-sections  of  flows  in  similar  hoppers  (see  text)  of  various  shapes; 
econe  =  30.680,  el/2  =  35.01°,  0PVramid  =  40°,  S  =  40°,  fiw  =  tan(20°),  N  =  U  =  20. 


dition,  tend  to  predict  flows  that  are  more  sensitive  to  wall  friction,  i.e.  the  material 
does  not  flow  as  well  near  the  wall,  [5],  [13],  [14]  ?  Unfortunately,  how  to  generalize 
Tresca’s  yield  condition  to  non-axisymmetric  flows  is  not  entirely  straightforward,  as 
one  looses  the  type  of  symmetry  properties  that  are  precisely  used  to  close  the  system 

2It  is  also  interesting  to  note  that  in  that  case,  the  steady  state  equations  are  always  hyperbolic. 
In  depth  numerical  calculations  under  a  Tresca  yield  condition  can  be  found  in  [4],  [13]. 
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in  the  usual  case;  see  [17],  p.29,  for  further  comments.  Accurate  flow  measurements 
(stress  and  velocity)  do  not  seem  to  be  available  at  this  time.  The  many  practical 
difficulties  related  to  precise  experiments  have  in  fact  been  at  the  origin  of  much  of 
the  analytical  work  in  this  field. 

6.  Conclusions.  We  have  shown  how  to  generalize  to  non-axisymmetric  con¬ 
tainers  the  notion  of  similarity  solutions  introduced  by  Jenike  [8]  in  the  case  of  purely 
conical  (or  wedge)  containers.  This  generalization  comes  at  the  price  of  having  to  solve 
a  nonlinear  Helmholtz  equation  on  a  part  of  the  sphere,  as  opposed  to  a  boundary 
problem  for  a  simple  system  of  ODEs  in  the  previous  cases. 

The  present  approach  applies  to  general  “conical”  domains,  in  the  mathematical 
sense,  i.e.,  any  domain  invariant  under  the  transformation  ( r,6,4 >)  i->  ( cr,8,4 >),  where 
c  >  0.  It  is  acknowledged  that,  clearly,  not  all  industrial  hoppers  satisfy  this  property. 
However,  failing  this,  no  similarity  solutions  are  expected  to  exist,  and  the  full  three- 
dimensional  equilibrium  equations  (2.1)  would  have  to  be  solved. 

Acknowledgments.  The  author  thanks  Tim  Kelley,  Matt  Matthews,  Tony  Royal, 
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